Antimicrobial resistance is becoming a worldwide problem
Many pathogens are becoming resistant for many different types of antibiotics
The institute of advanced studies (IAS) is doing research to the spreading of gonorrhoea in Amsterdam
For gonorrhoeae there might be a new antibiotic available within years
It is thus important to find a suitable management plan to prevent new resistance
Image(filename= "math_model.png")
plt.figure()
# plt.plot(t, sol_msm[1,:, 0], label='$S_{low}$')
# plt.plot(t, sol_msm[1,:, 1], label='$S_{high}$')
plt.plot(t, sol_hmw[1,:, 2], label='$I_{res, low}$')
plt.plot(t, sol_hmw[1,:, 3], label='$I_{res, high}$')
plt.plot(t, sol_hmw[1,:, 4], label='$I_{sen, low}$')
plt.plot(t, sol_hmw[1,:, 5], label='$I_{sen, high}$')
plt.legend(loc='best')
plt.xlabel('t (years)')
plt.title("Epidemic curve HMW group")
plt.ylabel("Fraction of population")
# plt.ylim(0, 1)
plt.xlim(0, 30)
plt.grid()
plt.show()
plt.figure()
ratio_sol_msm = (sol_msm[:,:, 2] + sol_msm[:,:, 3])/(sol_msm[:,:, 2] + sol_msm[:,:, 3] + sol_msm[:,:, 4] + sol_msm[:,:, 5])
plt.plot(t, ratio_sol_msm[1,:], label = "MSM", c = "g",linewidth=1)
plt.fill_between(t,ratio_sol_msm[0,:], ratio_sol_msm[2,:], alpha=0.5, facecolor = "g")
ratio_sol_hmw = (sol_hmw[:,:, 2] + sol_hmw[:,:, 3])/(sol_hmw[:,:, 2] + sol_hmw[:,:, 3] + sol_hmw[:,:, 4] + sol_hmw[:,:, 5])
plt.plot(t, ratio_sol_hmw[1,:], label = "HMW", c = "b",linewidth=1)
plt.fill_between(t,ratio_sol_hmw[0,:], ratio_sol_hmw[2,:], alpha=0.5, facecolor = "b")
plt.legend(loc='best')
plt.xlabel('t (years)')
plt.ylabel("Fraction of infected individuals")
plt.grid()
plt.ylim(0,1)
plt.show()
The parameters
# # create power law graph
# exp = 2.4
# seed = 123456789
# num = 1000
# graph = power_law_graph(exp, num, seed)
# or, create random network
degree = 10
num = 500
p = degree/num
seed = 1234
graph = random_network(p, num, seed)
# compute average degree
degrees = dict(graph.degree())
sum_of_edges = sum(degrees.values())
average_deg = sum_of_edges / len(degrees)
print(average_deg)
# compute r0 (TODO)
9.752
plt.figure()
pos = nx.spring_layout(graph)
nx.draw(graph, pos, node_color='b', node_size=10, with_labels=False)
plt.show()
# set person class to each node
for i in range(len(graph)):
graph.node[i]["Data"] = Individual(i)
# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)
# time steps
t = 0
t_tot = 20
time_step = 365
dt = 1 / (time_step)
# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]
# keep track of infecteds and steps
num_steps = int(t_tot / dt)
num_infected = np.zeros((2, num_steps))
num_res = np.zeros((2, num_steps))
p_inf = 0.8
# iterate over sextypes (MSM and HMW)
for sextype in range(2):
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
num_res[sextype,:], num_infected[sextype,:], graph = network_model(beta, tau, nu, mu, init, num_steps, copy.deepcopy(graph), True)
100%|██████████| 7300/7300 [00:44<00:00, 162.88it/s] 100%|██████████| 7300/7300 [00:42<00:00, 169.92it/s]
# plot number of infecteds
t = np.linspace(0,t_tot, num_steps)
plt.figure()
print("Mean number of infecteds:", np.mean(num_infected))
plt.plot(t, num_infected[0,:], label='all infecteds MSM')
plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
plt.plot(t, num_infected[1,:], label='all infecteds HMW')
plt.xlabel("Time (years)")
plt.ylabel("number of people")
plt.plot(t, num_res[1,:], label='resistant infecteds HMW')
plt.legend()
plt.show()
Mean number of infecteds: 288.4821917808219
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-280-cb51f78d178f> in <module>() 3 plt.figure() 4 print("Mean number of infecteds:", np.mean(num_infected)) ----> 5 plt.plot(t, num_infected[0,:], label='all infecteds MSM') 6 plt.plot(t, num_res[0,:], label='resistant infecteds MSM') 7 plt.plot(t, num_infected[1,:], label='all infecteds HMW') IndexError: too many indices for array
plt.figure()
plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
# plt.plot(num_res)
plt.xlabel("Time (years)")
plt.ylabel("Fraction of infected individuals")
plt.legend()
plt.show()
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) <ipython-input-281-9a4093e5464e> in <module>() 1 plt.figure() ----> 2 plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM') 3 plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW') 4 # plt.plot(num_res) 5 plt.xlabel("Time (years)") IndexError: too many indices for array
# set person class to each node
for i in range(len(graph)):
graph.node[i]["Data"] = Individual(i)
# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)
# time
t = 0
t_tot = 5
time_step = 365
dt = 1 / (time_step)
# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]
# choose MSM or HMW to model
sextype = 0
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
# initial state
_, _, graph = network_model(beta, tau, nu, mu, init, 0, graph, doInit=True, disable_progress=True)
# get values of disease status
values = []
for j in range(len(graph.nodes())):
if graph.node[j]["Data"].disease_status == 0:
values.append("Black")
elif graph.node[j]["Data"].disease_status == 1:
values.append("Blue")
elif graph.node[j]["Data"].disease_status == 2:
values.append("Red")
# create figure and set axes
fig = plt.figure(figsize = (10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1, = ax2.plot([], [], c = 'b', label='all infecteds MSM')
im2, = ax2.plot([], [], c = 'r', label='resistant infecteds MSM')
ax2.legend()
ax2.set_xlim(0,t_tot)
ax2.set_ylim(0,num)
# initialize
pos = nx.spring_layout(graph)
plt.close()
ani_step = 365//5
frames = (t_tot*365)//ani_step
num_infected = np.zeros(ani_step*frames)
num_res = np.zeros(ani_step*frames)
# update function for animation
def update(i):
if i<1:
temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, 1, graph, doInit=False, disable_progress=True)
else:
temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, ani_step, graph, doInit=False, disable_progress=True)
num_infected[i*ani_step:(i+1)*ani_step] = temp_infected
num_res[i*ani_step:(i+1)*ani_step] = temp_res
t = np.linspace(0,i*ani_step/365, i*ani_step)
im1.set_data(t, num_infected[:i*ani_step])
im2.set_data(t, num_res[:i*ani_step])
ax1.clear()
# solution at step i
values = []
for j in range(len(graph.nodes())):
if graph.node[j]["Data"].disease_status == 0:
values.append("Black")
elif graph.node[j]["Data"].disease_status == 1:
values.append("Blue")
elif graph.node[j]["Data"].disease_status == 2:
values.append("Red")
nx.draw(graph, pos=pos, node_color=values, node_size=10, ax=ax1, clim = 2)
# Scale plot ax
ax1.set_title("Frame %d: "%(i+1), fontweight="bold")
ax1.set_xticks([])
ax1.set_yticks([])
# animate!
animation.FuncAnimation(fig, update, frames=frames, interval=100, repeat=True)
25
The stochastic network model is based on an existing model by Kretzschmar, M. et al. (1996) [4]
The system consists of Persons defined by:
During each timestep, Partnerships are formed between two Persons. A Partnership is defined by:
The sexual contact network is therefore dynamic in time
Each timestep of the model consists of seven parts:
The resistant disease enters the system through treatment; treatment will result in contracting the resistant disease with probability 0.01%
fig
fig
fig
fig
Although the disease might seem quite infectious, simulations have shown $R_0=0.90 \pm 2.93$ (N=1,000, 95% CI). The spread of the disease is very dependent on where the disease is started. If a person with high sexual activity is infected, it will quickly spread throughout the entire core group whereas a person with low sexual activity may only affect a select few different partners
anim
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-211-027b5d0310db> in <module>() ----> 1 anim /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in __call__(self, result) 255 self.start_displayhook() 256 self.write_output_prompt() --> 257 format_dict, md_dict = self.compute_format_data(result) 258 self.update_user_ns(result) 259 self.fill_exec_result(result) /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in compute_format_data(self, result) 149 150 """ --> 151 return self.shell.display_formatter.format(result) 152 153 # This can be set to True by the write_output_prompt method in a subclass /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude) 178 md = None 179 try: --> 180 data = formatter(obj) 181 except: 182 # FIXME: log the exception <decorator-gen-9> in __call__(self, obj) /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs) 222 """show traceback on failed format call""" 223 try: --> 224 r = method(self, *args, **kwargs) 225 except NotImplementedError: 226 # don't warn on NotImplementedErrors /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in __call__(self, obj) 343 method = get_real_method(obj, self.print_method) 344 if method is not None: --> 345 return method() 346 return None 347 else: /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _repr_html_(self) 1484 return self.to_html5_video() 1485 elif fmt == 'jshtml': -> 1486 return self.to_jshtml() 1487 1488 /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in to_jshtml(self, fps, embed_frames, default_mode) 1467 self.save(f.name, writer=HTMLWriter(fps=fps, 1468 embed_frames=embed_frames, -> 1469 default_mode=default_mode)) 1470 # Re-open and get content 1471 with open(f.name) as fobj: /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs) 1257 for anim, d in zip(all_anim, data): 1258 # TODO: See if turning off blit is really necessary -> 1259 anim._draw_next_frame(d, blit=False) 1260 writer.grab_frame(**savefig_kwargs) 1261 /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_next_frame(self, framedata, blit) 1294 # post- draw, as well as the drawing of the frame itself. 1295 self._pre_draw(framedata, blit) -> 1296 self._draw_frame(framedata) 1297 self._post_draw(framedata, blit) 1298 /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_frame(self, framedata) 1812 # Call the func with framedata and args. If blitting is desired, 1813 # func needs to return a sequence of any artists that were modified. -> 1814 self._drawn_artists = self._func(framedata, *self._args) 1815 if self._blit: 1816 if self._drawn_artists is None: <ipython-input-210-e04a745e332f> in plot_graph(i) 29 if i > 0: 30 for j in range(365): ---> 31 s.time_step() 32 33 persons = set() ~/Documents/GitHub/CSS2018/Code/stochastic_network_model.py in time_step(self) 645 if np.random.random() < self.recovery_rate_asymptomatic_men: 646 p.cure(self) --> 647 elif p.gender == 1 and p.disease_status in [1,3]: 648 if np.random.random() < self.recovery_rate_symptomatic_women: 649 p.cure(self) KeyboardInterrupt:
[1] Fingerhuth, S. M., Bonhoeffer, S., Low, N., & Althaus, C. L. (2016). Antibiotic-resistant Neisseria gonorrhoeae spread faster with more treatment, not more sexual partners. PLoS pathogens, 12(5), e1005611.
[2] De, P., Singh, A. E., Wong, T., Yacoub, W., & Jolly, A. M. (2004). Sexual network analysis of a gonorrhoea outbreak. Sexually transmitted infections, 80(4), 280-285.
[3] Bansal, S., Grenfell, B. T., & Meyers, L. A. (2007). When individual behaviour matters: homogeneous and network models in epidemiology. Journal of the Royal Society Interface, 4(16), 879-891.
[4] Kretzschmar, M., van Duynhoven, Y. T., & Severijnen, A. J. (1996). Modeling prevention strategies for gonorrhea and chlamydia using stochastic network simulations. American Journal of Epidemiology, 144(3), 306-317.